library(lme4) # pour la modélisation multiniveau
library(multilevelTools) # pour l'analyse des modèles
library(lmerTest) # pour l'évaluation des modèles
library(sf) # pour les objets spatiaux (cartes)
library(tidyverse) # pour la manipulation des données
library(readxl) # pour ouvrir les fichiers xlsx
library(ggrepel) # pour des étiquettes de graphiques propres
Il s’agit des vignobles de Bourgogne. On a dans un tableau les caractéristiques de 2391 vignobles en termes de:
vignobles <- st_read("Data/Geo/Vignobles_prix.shp",
stringsAsFactors = T) %>%
st_transform(crs = "EPSG:2154")
## Reading layer `Vignobles_prix' from data source
## `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Geo/Vignobles_prix.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2391 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 826631.6 ymin: 6645988 xmax: 851465.7 ymax: 6690810
## Projected CRS: RGF93 v1 / Lambert-93
summary(vignobles)
## Concat LIBCOM
## 1AUXEY-DURESSESCORLOIN : 1 MEURSAULT : 158
## 1AUXEY-DURESSESDERRIERE LE BOIS: 1 GEVREY-CHAMBERTIN : 128
## 1AUXEY-DURESSESEN CHAIGNUT : 1 BEAUNE : 126
## 1AUXEY-DURESSESEN FRIVOLE : 1 NUITS-SAINT-GEORGES: 121
## 1AUXEY-DURESSESEN MORVIN : 1 SAINT-AUBIN : 121
## 1AUXEY-DURESSESEN POCHENOT : 1 (Other) :1727
## (Other) :2385 NA's : 10
## NOM NIVEAU SURFACE SCORE_b
## LE VILLAGE : 26 Bourgogne : 624 Min. :0.0000 Min. :18.40
## LES CRAIS : 10 Coteaux b. : 319 1st Qu.:0.1300 1st Qu.:67.90
## LES CRAS : 8 Grand cru : 37 Median :0.3100 Median :73.20
## SOUS LA VELLE: 8 Premier cru: 410 Mean :0.4773 Mean :71.89
## LES PERRIERES: 7 Village :1001 3rd Qu.:0.6100 3rd Qu.:76.60
## (Other) :2331 Max. :7.2000 Max. :94.22
## NA's : 1
## Source PrixMoy geometry
## guide_hachette: 416 Min. : 5.0 MULTIPOLYGON :2391
## interpolation :1942 1st Qu.: 23.0 epsg:2154 : 0
## wine_search : 33 Median : 36.0 +proj=lcc ...: 0
## Mean : 154.4
## 3rd Qu.: 87.0
## Max. :25960.0
## NA's :4
vignobles %>%
dplyr::filter(LIBCOM == "FLAGEY-ECHEZEAUX" & NIVEAU == "Premier cru")
Il s’agit des communes de Côte d’Or. On a dans un tableau les caractéristiques de 31 communes en termes de:
communes <- st_read("Data/Geo/communes.shp",
stringsAsFactors = T) %>%
st_transform(crs = "EPSG:2154")
## Reading layer `Communes' from data source
## `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Geo/Communes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 31 features and 20 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 826428 ymin: 6645899 xmax: 855115 ymax: 6691133
## Projected CRS: RGF93 v1 / Lambert-93
summary(communes)
## gid id_geofla code_comm insee_com
## Min. : 563 Min. : 563 010 : 1 21010 : 1
## 1st Qu.: 6744 1st Qu.: 6744 037 : 1 21037 : 1
## Median :14409 Median :14409 054 : 1 21054 : 1
## Mean :11410 Mean :11410 110 : 1 21110 : 1
## 3rd Qu.:14602 3rd Qu.:14602 133 : 1 21133 : 1
## Max. :14749 Max. :14749 150 : 1 21150 : 1
## (Other):25 (Other):25
## nom_comm statut x_chf_lieu y_chf_lieu
## ALOXE-CORTON : 1 Chef-lieu canton: 3 Min. :8291 Min. :66472
## AUXEY-DURESSES : 1 Commune simple :27 1st Qu.:8349 1st Qu.:66570
## BEAUNE : 1 Sous-préfecture : 1 Median :8431 Median :66660
## BROCHON : 1 Mean :8419 Mean :66676
## CHAMBOLLE-MUSIGNY : 1 3rd Qu.:8486 3rd Qu.:66774
## CHASSAGNE-MONTRACHET: 1 Max. :8516 Max. :66897
## (Other) :25
## x_centroid y_centroid z_moyen superficie population
## Min. :8288 Min. :66477 Min. :222.0 Min. : 90 Min. : 0.20
## 1st Qu.:8350 1st Qu.:66576 1st Qu.:254.0 1st Qu.: 687 1st Qu.: 0.35
## Median :8450 Median :66661 Median :293.0 Median : 891 Median : 0.60
## Mean :8418 Mean :66677 Mean :303.6 Mean :1147 Mean : 2.21
## 3rd Qu.:8482 3rd Qu.:66774 3rd Qu.:349.0 3rd Qu.:1267 3rd Qu.: 1.30
## Max. :8514 Max. :66899 Max. :464.0 Max. :3619 Max. :21.90
##
## code_cant code_arr code_dept nom_dept code_reg nom_region
## 05 :8 1:23 21:31 COTE-D'OR:31 26:31 BOURGOGNE:31
## 24 :7 2: 8
## 15 :6
## 23 :5
## 06 :2
## 38 :1
## (Other):2
## Cote geometry
## Côte D'or : 3 POLYGON :31
## Côte de Beane:16 epsg:2154 : 0
## Côte de Nuit :12 +proj=lcc ...: 0
##
##
##
##
Avec ggplot, utilisation de geom_sf pour
les objets spatiaux.
ggplot(data = vignobles) +
geom_sf(aes(fill = log(PrixMoy)), lwd=0.01, colour="white") +
scale_fill_viridis_c(option = "viridis")
Pour des étiquettes réactives, on utilise
geom_text_repel:
ggplot(data = communes) +
geom_sf(aes(fill = Cote), lwd=0.01, colour="white") +
geom_text_repel(data = communes[communes$population > 1,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
Info sur la distribution des prix, des surfaces, des qualités physiques des vignobles, etc.
Vignobles_par_communes <- as_tibble(vignobles) %>%
group_by(LIBCOM) %>%
summarise(aire_vignoble = sum(SURFACE),
MoyPond_b = sum(SCORE_b * SURFACE) / sum(SURFACE),
N_Parcelles = n(),
Prix_moyen = mean(PrixMoy,na.rm=T),
Prix_median = median(PrixMoy,na.rm=T),
)
Info sur les appellations d’origine contrôlée:
Vignobles_par_communes_par_AOC <- as_tibble(vignobles) %>%
group_by(LIBCOM, NIVEAU) %>%
summarise(N_Parcelles = n()
) %>%
pivot_wider(names_from = NIVEAU, values_from = N_Parcelles) %>%
replace_na(list(`Bourgogne` = 0,
`Grand cru`=0,
`Premier cru`=0,
`Village`=0,
`Coteaux b.`=0))
## `summarise()` has grouped output by 'LIBCOM'. You can override using the
## `.groups` argument.
communes_final <- left_join(Vignobles_par_communes,
Vignobles_par_communes_par_AOC, by="LIBCOM") %>%
left_join(communes,., by=c("nom_comm"="LIBCOM")) %>%
mutate(share_Bourgogne = `Bourgogne` / N_Parcelles * 100,
share_GdCru = `Grand cru` / N_Parcelles * 100,
share_PmCru = `Premier cru` / N_Parcelles * 100,
share_Village = `Village` / N_Parcelles * 100,
share_Coteau = `Coteaux b.` / N_Parcelles * 100,
share_PmGrCru = (`Premier cru` + `Grand cru`) / N_Parcelles * 100,
share_PmGrCruVill = (`Premier cru` + `Grand cru` + `Village`) / N_Parcelles * 100)
head(communes_final)
# Prix moyen des vins estimés par commune
ggplot(data = communes_final) +
geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white") +
scale_fill_gradient(low = "white", high="red") +
geom_text_repel(data = communes_final[communes_final$Prix_moyen > 100 ,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
# Qualité moyenne des parcelles par commune
ggplot(data = communes_final) +
geom_sf(aes(fill = MoyPond_b), lwd=0.01, colour="white") +
scale_fill_gradient(low = "white", high="goldenrod4") +
geom_text_repel(data = communes_final[communes_final$MoyPond_b > 80
| communes_final$MoyPond_b < 65 ,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
# Proportion de parcelles en grand cru par commune
ggplot(data = communes_final) +
geom_sf(aes(fill = share_GdCru), lwd=0.01, colour="white") +
scale_fill_gradient(low = "white", high="goldenrod") +
geom_text_repel(data = communes_final[communes_final$share_GdCru > 8 ,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
# Proportion de parcelles en premier ou grand cru par commune
ggplot(data = communes_final) +
geom_sf(aes(fill = share_PmGrCru), lwd=0.01, colour="white") +
scale_fill_gradient(low = "white", high="gold") +
geom_text_repel(data = communes_final[communes_final$share_PmGrCru > 30 ,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
communes_final_tib <- as_tibble(communes_final) %>%
select(-geometry)
vignobles_final <- left_join(vignobles,communes_final_tib,
by=c("LIBCOM"= "nom_comm"))
summary(vignobles_final)
## Concat LIBCOM
## 1AUXEY-DURESSESCORLOIN : 1 MEURSAULT : 158
## 1AUXEY-DURESSESDERRIERE LE BOIS: 1 GEVREY-CHAMBERTIN : 128
## 1AUXEY-DURESSESEN CHAIGNUT : 1 BEAUNE : 126
## 1AUXEY-DURESSESEN FRIVOLE : 1 NUITS-SAINT-GEORGES: 121
## 1AUXEY-DURESSESEN MORVIN : 1 SAINT-AUBIN : 121
## 1AUXEY-DURESSESEN POCHENOT : 1 (Other) :1727
## (Other) :2385 NA's : 10
## NOM NIVEAU SURFACE SCORE_b
## LE VILLAGE : 26 Bourgogne : 624 Min. :0.0000 Min. :18.40
## LES CRAIS : 10 Coteaux b. : 319 1st Qu.:0.1300 1st Qu.:67.90
## LES CRAS : 8 Grand cru : 37 Median :0.3100 Median :73.20
## SOUS LA VELLE: 8 Premier cru: 410 Mean :0.4773 Mean :71.89
## LES PERRIERES: 7 Village :1001 3rd Qu.:0.6100 3rd Qu.:76.60
## (Other) :2331 Max. :7.2000 Max. :94.22
## NA's : 1
## Source PrixMoy gid id_geofla
## guide_hachette: 416 Min. : 5.0 Min. : 563 Min. : 563
## interpolation :1942 1st Qu.: 23.0 1st Qu.: 6720 1st Qu.: 6720
## wine_search : 33 Median : 36.0 Median :14392 Median :14392
## Mean : 154.4 Mean :10484 Mean :10484
## 3rd Qu.: 87.0 3rd Qu.:14597 3rd Qu.:14597
## Max. :25960.0 Max. :14749 Max. :14749
## NA's :4 NA's :10 NA's :10
## code_comm insee_com statut x_chf_lieu
## 412 : 158 21412 : 158 Chef-lieu canton: 257 Min. :8291
## 295 : 128 21295 : 128 Commune simple :1998 1st Qu.:8342
## 054 : 126 21054 : 126 Sous-préfecture : 126 Median :8396
## 464 : 121 21464 : 121 NA's : 10 Mean :8405
## 541 : 121 21541 : 121 3rd Qu.:8480
## (Other):1727 (Other):1727 Max. :8516
## NA's : 10 NA's : 10 NA's :10
## y_chf_lieu x_centroid y_centroid z_moyen superficie
## Min. :66472 Min. :8288 Min. :66477 Min. :222.0 Min. : 90
## 1st Qu.:66555 1st Qu.:8338 1st Qu.:66553 1st Qu.:256.0 1st Qu.: 765
## Median :66641 Median :8395 Median :66642 Median :300.0 Median :1006
## Mean :66654 Mean :8404 Mean :66655 Mean :306.2 Mean :1353
## 3rd Qu.:66758 3rd Qu.:8477 3rd Qu.:66751 3rd Qu.:355.0 3rd Qu.:1962
## Max. :66897 Max. :8514 Max. :66899 Max. :464.0 Max. :3619
## NA's :10 NA's :10 NA's :10 NA's :10 NA's :10
## population code_cant code_arr code_dept nom_dept
## Min. : 0.200 05 :688 1 :1795 21 :2381 COTE-D'OR:2381
## 1st Qu.: 0.400 15 :515 2 : 586 NA's: 10 NA's : 10
## Median : 0.700 23 :457 NA's: 10
## Mean : 2.385 24 :374
## 3rd Qu.: 1.600 06 :150
## Max. :21.900 (Other):197
## NA's :10 NA's : 10
## code_reg nom_region Cote aire_vignoble
## 26 :2381 BOURGOGNE:2381 Côte D'or : 152 Min. : 5.25
## NA's: 10 NA's : 10 Côte de Beane:1421 1st Qu.: 23.39
## Côte de Nuit : 808 Median : 44.57
## NA's : 10 Mean : 45.51
## 3rd Qu.: 58.18
## Max. :101.94
## NA's :10
## MoyPond_b N_Parcelles Prix_moyen Prix_median
## Min. :58.19 Min. : 8.00 Min. : 16.55 Min. : 17.0
## 1st Qu.:69.98 1st Qu.: 63.00 1st Qu.: 26.11 1st Qu.: 26.0
## Median :71.39 Median : 88.00 Median : 35.48 Median : 34.5
## Mean :71.07 Mean : 93.89 Mean : 154.16 Mean : 112.9
## 3rd Qu.:73.59 3rd Qu.:121.00 3rd Qu.: 140.79 3rd Qu.: 84.0
## Max. :87.55 Max. :158.00 Max. :2151.34 Max. :1447.0
## NA's :10 NA's :10 NA's :10 NA's :10
## Bourgogne Grand cru Premier cru Village
## Min. : 2.00 Min. :0.000 Min. : 0.00 Min. : 2.00
## 1st Qu.:12.00 1st Qu.:0.000 1st Qu.: 8.00 1st Qu.:27.00
## Median :21.00 Median :0.000 Median :17.00 Median :39.00
## Mean :24.49 Mean :1.351 Mean :16.81 Mean :38.15
## 3rd Qu.:34.00 3rd Qu.:2.000 3rd Qu.:26.00 3rd Qu.:48.00
## Max. :55.00 Max. :9.000 Max. :42.00 Max. :62.00
## NA's :10 NA's :10 NA's :10 NA's :10
## Coteaux b. share_Bourgogne share_GdCru share_PmCru
## Min. : 0.00 Min. :11.11 Min. : 0.000 Min. : 0.000
## 1st Qu.: 5.00 1st Qu.:19.05 1st Qu.: 0.000 1st Qu.: 7.207
## Median :10.00 Median :26.98 Median : 0.000 Median :19.672
## Mean :13.09 Mean :26.08 Mean : 1.554 Mean :17.136
## 3rd Qu.:16.00 3rd Qu.:31.40 3rd Qu.: 2.500 3rd Qu.:24.324
## Max. :41.00 Max. :75.00 Max. :11.111 Max. :44.444
## NA's :10 NA's :10 NA's :10 NA's :10
## share_Village share_Coteau share_PmGrCru share_PmGrCruVill
## Min. :11.54 Min. : 0.000 Min. : 0.000 Min. :25.00
## 1st Qu.:32.77 1st Qu.: 4.688 1st Qu.: 7.207 1st Qu.:46.84
## Median :39.24 Median :11.570 Median :23.529 Median :63.22
## Mean :41.83 Mean :13.398 Mean :18.690 Mean :60.52
## 3rd Qu.:48.98 3rd Qu.:17.857 3rd Qu.:27.473 3rd Qu.:73.75
## Max. :76.54 Max. :44.444 Max. :55.556 Max. :86.49
## NA's :10 NA's :10 NA's :10 NA's :10
## geometry
## MULTIPOLYGON :2391
## epsg:2154 : 0
## +proj=lcc ...: 0
##
##
##
##
summary(vignobles_final$PrixMoy)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 5.0 23.0 36.0 154.4 87.0 25960.0 4
ggplot(data=vignobles_final, aes(x=PrixMoy)) +
geom_histogram(fill="coral1")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).
# Passage en log pour une distribution plus normale
vignobles_final$LogPrix <- log(vignobles_final$PrixMoy)
ggplot(data=vignobles_final, aes(x=LogPrix)) +
geom_histogram(fill="coral3") +
geom_vline(aes(xintercept=mean(LogPrix)),
color="black", linetype="dashed", size=1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 2391 rows containing missing values (`geom_vline()`).
#Score qualité geographique (pente, geologie, pedologie, precipitations):
ggplot(data=vignobles_final, aes(x=SCORE_b)) +
geom_histogram(fill="brown1")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Surface:
ggplot(data=vignobles_final, aes(x=SURFACE)) +
geom_histogram(fill="brown2")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Niveau d'aoc:
positions <- c("Coteaux b.", "Bourgogne", "Village",
"Premier cru", "Grand cru")
ggplot(data=vignobles_final, aes(x=NIVEAU)) +
geom_bar(fill="brown3") +
scale_x_discrete(limits = positions)
#Source info prix:
ggplot(data=vignobles_final, aes(x=Source)) +
geom_bar(fill="brown4")
#Population de la commune
ggplot(data=communes_final, aes(x=population)) +
geom_histogram(fill="aquamarine1", bins = 5)
#Part de grand et premier cru de la commune
ggplot(data=communes_final, aes(x=share_PmGrCru)) +
geom_histogram(fill="aquamarine3", bins = 5)
#Cote de nuit ou cote de beaune:
pos_cote <- c("Côte de Nuit", "Côte de Beane", "Côte D'or")
ggplot(data=communes_final, aes(x=Cote)) +
geom_bar(fill="aquamarine4") +
scale_x_discrete(limits = pos_cote)
#Qualite moyenne de la commune
ggplot(data=communes_final, aes(x=MoyPond_b)) +
geom_histogram(fill="darkslategrey", bins = 5)
vignobles_final_cr <- vignobles_final %>%
mutate(LogPrix = scale(as.numeric(LogPrix)),
Surface = scale(as.numeric(SURFACE)),
Quality = scale(as.numeric(SCORE_b)),
AOC = fct_reorder(NIVEAU, SCORE_b),
Cote = fct_reorder(Cote, SCORE_b),
MeanQuality = scale(as.numeric(MoyPond_b)),
ShareGdCru = scale(as.numeric(share_GdCru))) %>%
select(LogPrix, Surface, Quality, AOC,
MeanQuality, ShareGdCru,
LIBCOM, Source, Cote, Concat) %>%
na.omit
summary_stat <- data.frame()
L1 <- lm(LogPrix~Quality+Surface+AOC+Source, data=vignobles_final_cr)
summary(L1)
##
## Call:
## lm(formula = LogPrix ~ Quality + Surface + AOC + Source, data = vignobles_final_cr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0696 -0.6275 -0.2024 0.4016 3.9311
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.71893 0.09121 -7.882 4.88e-15 ***
## Quality -0.01648 0.03248 -0.507 0.611999
## Surface -0.07580 0.01951 -3.885 0.000105 ***
## AOCBourgogne -0.15699 0.07522 -2.087 0.036995 *
## AOCVillage -0.20275 0.08539 -2.374 0.017655 *
## AOCPremier cru -0.05133 0.11279 -0.455 0.649102
## AOCGrand cru 2.33775 0.46750 5.001 6.14e-07 ***
## Sourceinterpolation 0.98998 0.05650 17.520 < 2e-16 ***
## Sourcewine_search 0.96382 0.46841 2.058 0.039735 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8659 on 2368 degrees of freedom
## Multiple R-squared: 0.2521, Adjusted R-squared: 0.2496
## F-statistic: 99.8 on 8 and 2368 DF, p-value: < 2.2e-16
mnnull <- lmer(LogPrix ~ 1 + (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mnnull)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: LogPrix ~ 1 + (1 | LIBCOM)
## Data: vignobles_final_cr
##
## AIC BIC logLik deviance df.resid
## 4487.1 4504.4 -2240.5 4481.1 2374
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4605 -0.3560 0.0610 0.4361 4.9734
##
## Random effects:
## Groups Name Variance Std.Dev.
## LIBCOM (Intercept) 0.8638 0.9294
## Residual 0.3611 0.6009
## Number of obs: 2377, groups: LIBCOM, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.1520 0.1677 30.9052 0.907 0.372
vrn2 <- summary(mnnull)$varcor[1][[1]][[1]]
vrn1 <- as.data.frame(summary(mnnull)$varcor)[2,]$vcov
icc <- vrn2 / (vrn1 + vrn2)
icc
## [1] 0.7051989
Ou plus simplement:
performance::icc(mnnull)[[1]]
## [1] 0.7051989
Dans ce modèle vide à 2 niveaux, 29.5% de la variance est prise en charge par l’appartenance au niveau 2.
# Pour résumer:
summary_stat["MNNull","N"] <- dim(mnnull@frame)[1]
summary_stat["MNNull","LogLik" ] <- summary(mnnull)$AIC[3][[1]]
summary_stat["MNNull", "AIC" ] <-summary(mnnull)$AIC[1][[1]]
summary_stat["MNNull", "BIC" ] <- summary(mnnull)$AIC[2][[1]]
summary_stat["MNNull", "Deviance" ] <- summary(mnnull)$AIC[4][[1]]
summary_stat["MNNull", "Total_Var" ] <- vrn2 + vrn1
summary_stat["MNNull", "ICC" ] <- performance::icc(mnnull)[[1]]
summary_stat["MNNull", "InterClass" ] <- vrn1 / (vrn2 + vrn1)
summary_stat
Extraction des résidus avec ranef, puis jointure et
cartographie.
ranef <- ranef(mnnull)$LIBCOM %>%
mutate(LIBCOM = as.character(rownames(.)),
residuCom = as.numeric(`(Intercept)`))
mnnull_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))
ggplot(data = mnnull_ranef) +
geom_sf(aes(fill = residuCom), lwd=0.01, colour="white") +
scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
geom_text_repel(data = mnnull_ranef[abs(mnnull_ranef$residuCom) > 1,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
À comparer avec la carte originale des prix moyens par commune = différence de centrage-réduction.
ggplot(data = communes_final) +
geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white") +
scale_fill_gradient(low = "white", high="red") +
geom_text_repel(data = communes_final[communes_final$Prix_moyen > 500 ,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
mn1 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
(1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: LogPrix ~ Quality + Surface + Source + AOC + (1 | LIBCOM)
## Data: vignobles_final_cr
##
## AIC BIC logLik deviance df.resid
## 3038.4 3101.9 -1508.2 3016.4 2366
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4448 -0.5237 -0.0687 0.4612 4.6352
##
## Random effects:
## Groups Name Variance Std.Dev.
## LIBCOM (Intercept) 0.7839 0.8854
## Residual 0.1937 0.4401
## Number of obs: 2377, groups: LIBCOM, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.50786 0.16599 36.31713 -3.060 0.004149 **
## Quality -0.04720 0.01704 2347.41662 -2.770 0.005644 **
## Surface 0.04862 0.01028 2347.53697 4.728 2.40e-06 ***
## Sourceinterpolation 0.89024 0.02919 2346.53970 30.498 < 2e-16 ***
## Sourcewine_search 1.90564 0.24399 2347.88348 7.810 8.53e-15 ***
## AOCBourgogne -0.07872 0.03881 2346.73517 -2.028 0.042646 *
## AOCVillage -0.17123 0.04409 2346.68467 -3.884 0.000106 ***
## AOCPremier cru -0.07690 0.05793 2346.38460 -1.327 0.184518
## AOCGrand cru 0.46195 0.24339 2347.85577 1.898 0.057822 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Qualty Surfac Srcntr Srcwn_ AOCBrg AOCVll AOCPrc
## Quality 0.152
## Surface -0.010 0.229
## Sorcntrpltn -0.168 0.072 0.212
## Sorcwn_srch -0.006 0.060 0.191 0.094
## AOCBourgogn -0.197 -0.584 -0.073 -0.028 -0.028
## AOCVillage -0.229 -0.720 -0.086 0.087 -0.022 0.817
## AOCPremircr -0.237 -0.745 -0.096 0.234 -0.014 0.743 0.848
## AOCGrandcru -0.067 -0.292 -0.258 0.002 -0.908 0.232 0.268 0.271
summary_stat["MNVar1Niv","N"] <- dim(mn1@frame)[1]
summary_stat["MNVar1Niv","LogLik" ] <- summary(mn1)$AIC[3][[1]]
summary_stat["MNVar1Niv", "AIC" ] <-summary(mn1)$AIC[1][[1]]
summary_stat["MNVar1Niv", "BIC" ] <- summary(mn1)$AIC[2][[1]]
summary_stat["MNVar1Niv", "Deviance" ] <- summary(mn1)$AIC[4][[1]]
summary_stat["MNVar1Niv", "Total_Var" ] <-
summary(mn1)$varcor[1][[1]][[1]] + as.data.frame(summary(mn1)$varcor)[2,]$vcov
summary_stat["MNVar1Niv", "ICC" ] <- performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "InterClass" ] <- 1 - performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "PseudoR2" ] <-
(1 - (summary_stat["MNVar1Niv", "Total_Var" ] / summary_stat["MNNull", "Total_Var" ] ))
summary_stat
Résidus de niveau 1
mn1_residuals <- cbind(vignobles_final_cr, residuals(mn1))
colnames(mn1_residuals)[length(colnames(vignobles_final_cr))] <- "residusN1"
summary(mn1_residuals$residusN1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.39613 -0.23049 -0.03025 0.00000 0.20294 2.03984
ggplot(data = mn1_residuals) +
geom_sf(aes(fill = residusN1), lwd=0.01, colour="white") +
scale_fill_gradient2(low = "darkgreen", mid = "white", high = "darkorange") +
theme_minimal()
Résidus de niveau 2
ranef <- ranef(mn1)$LIBCOM %>%
mutate(LIBCOM = as.character(rownames(.)),
residuCom = as.numeric(`(Intercept)`))
mn1_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))
ggplot(data = mn1_ranef) +
geom_sf(aes(fill = residuCom), lwd=0.01, colour="white") +
scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
geom_text_repel(data = mn1_ranef[abs(mn1_ranef$residuCom) > 1,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
mn2 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
MeanQuality + ShareGdCru + Cote +
(1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## LogPrix ~ Quality + Surface + Source + AOC + MeanQuality + ShareGdCru +
## Cote + (1 | LIBCOM)
## Data: vignobles_final_cr
##
## AIC BIC logLik deviance df.resid
## 2999.6 3086.2 -1484.8 2969.6 2362
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4420 -0.5218 -0.0649 0.4570 4.6376
##
## Random effects:
## Groups Name Variance Std.Dev.
## LIBCOM (Intercept) 0.1684 0.4104
## Residual 0.1937 0.4401
## Number of obs: 2377, groups: LIBCOM, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.05896 0.11393 42.28250 -9.295 8.91e-12 ***
## Quality -0.04788 0.01704 2346.30069 -2.809 0.005004 **
## Surface 0.04860 0.01028 2350.17563 4.727 2.41e-06 ***
## Sourceinterpolation 0.89129 0.02919 2347.54433 30.537 < 2e-16 ***
## Sourcewine_search 1.88415 0.24384 2353.35431 7.727 1.62e-14 ***
## AOCBourgogne -0.07774 0.03880 2347.94052 -2.003 0.045241 *
## AOCVillage -0.17095 0.04408 2347.70878 -3.878 0.000108 ***
## AOCPremier cru -0.07587 0.05793 2347.03637 -1.310 0.190419
## AOCGrand cru 0.47943 0.24325 2352.82532 1.971 0.048845 *
## MeanQuality -0.04275 0.08039 30.57464 -0.532 0.598741
## ShareGdCru 0.36251 0.08666 30.27123 4.183 0.000227 ***
## CoteCôte D'or 0.62786 0.26975 31.72140 2.328 0.026477 *
## CoteCôte de Nuit 1.14890 0.17176 30.09307 6.689 2.04e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
ShareGdCru).summary_stat["MNVar2Niv","N"] <- dim(mn2@frame)[1]
summary_stat["MNVar2Niv","LogLik" ] <- summary(mn2)$AIC[3][[1]]
summary_stat["MNVar2Niv", "AIC" ] <-summary(mn2)$AIC[1][[1]]
summary_stat["MNVar2Niv", "BIC" ] <- summary(mn2)$AIC[2][[1]]
summary_stat["MNVar2Niv", "Deviance" ] <- summary(mn2)$AIC[4][[1]]
summary_stat["MNVar2Niv", "Total_Var" ] <-
summary(mn2)$varcor[1][[1]][[1]] + as.data.frame(summary(mn2)$varcor)[2,]$vcov
summary_stat["MNVar2Niv", "ICC" ] <- performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "InterClass" ] <- 1 - performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "PseudoR2" ] <-
(1 - (summary_stat["MNVar2Niv", "Total_Var" ] / summary_stat["MNNull", "Total_Var" ] ))
summary_stat
Résidus de niveau 1
mn2_residuals <- cbind(vignobles_final_cr, residuals(mn2))
colnames(mn2_residuals)[length(colnames(vignobles_final_cr))] <- "residusN1"
ggplot(data = mn2_residuals) +
geom_sf(aes(fill = residusN1), lwd=0.01, colour="white") +
scale_fill_gradient2(low = "darkgreen", mid = "white", high = "darkorange") +
theme_minimal()
Résidus de niveau 2
ranef <- ranef(mn2)$LIBCOM %>%
mutate(LIBCOM = as.character(rownames(.)),
residuCom = as.numeric(`(Intercept)`))
mn2_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))
ggplot(data = mn2_ranef) +
geom_sf(aes(fill = residuCom), lwd=0.01, colour="white") +
scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
geom_text_repel(data = mn2_ranef[abs(mn2_ranef$residuCom) > 1,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
ggplot(data = vignobles_final_cr,
aes(x = Quality,
y = LogPrix,
col = LIBCOM))+
geom_point(size = 1,
alpha = .7,
position = "jitter")+
geom_smooth(method = lm,
se = T,
size = 1.5,
linetype = 1,
alpha = .2) +
labs(title = "Relation entre qualité et prix selon l'appartenance communale") +
theme(legend.position="right")
## `geom_smooth()` using formula = 'y ~ x'
mn3 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
MeanQuality + ShareGdCru + Cote +
(Quality | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula:
## LogPrix ~ Quality + Surface + Source + AOC + MeanQuality + ShareGdCru +
## Cote + (Quality | LIBCOM)
## Data: vignobles_final_cr
##
## AIC BIC logLik deviance df.resid
## 2796.9 2895.0 -1381.4 2762.9 2360
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7412 -0.5349 -0.0108 0.4587 4.4033
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## LIBCOM (Intercept) 0.16457 0.4057
## Quality 0.03362 0.1834 0.23
## Residual 0.17257 0.4154
## Number of obs: 2377, groups: LIBCOM, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.167e+00 1.122e-01 4.315e+01 -10.403 2.46e-13 ***
## Quality -1.006e-01 4.025e-02 4.262e+01 -2.499 0.016381 *
## Surface 4.033e-02 9.926e-03 2.348e+03 4.064 4.99e-05 ***
## Sourceinterpolation 8.805e-01 2.797e-02 2.335e+03 31.483 < 2e-16 ***
## Sourcewine_search 2.145e+00 2.518e-01 2.233e+03 8.519 < 2e-16 ***
## AOCBourgogne -4.226e-02 3.796e-02 2.346e+03 -1.113 0.265788
## AOCVillage -1.083e-01 4.447e-02 2.346e+03 -2.435 0.014961 *
## AOCPremier cru 3.637e-03 5.951e-02 2.347e+03 0.061 0.951272
## AOCGrand cru 4.217e-01 2.521e-01 2.251e+03 1.673 0.094523 .
## MeanQuality -4.740e-02 7.855e-02 3.139e+01 -0.603 0.550536
## ShareGdCru 3.694e-01 8.405e-02 3.036e+01 4.395 0.000125 ***
## CoteCôte D'or 7.577e-01 2.616e-01 3.176e+01 2.896 0.006783 **
## CoteCôte de Nuit 1.298e+00 1.663e-01 3.004e+01 7.807 1.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
summary_stat["MNVar3Niv","N"] <- dim(mn3@frame)[1]
summary_stat["MNVar3Niv","LogLik" ] <- summary(mn3)$AIC[3][[1]]
summary_stat["MNVar3Niv", "AIC" ] <-summary(mn3)$AIC[1][[1]]
summary_stat["MNVar3Niv", "BIC" ] <- summary(mn3)$AIC[2][[1]]
summary_stat["MNVar3Niv", "Deviance" ] <- summary(mn3)$AIC[4][[1]]
summary_stat["MNVar3Niv", "Total_Var" ] <-
summary(mn3)$varcor[1][[1]][[1]] + as.data.frame(summary(mn3)$varcor)[2,]$vcov
summary_stat["MNVar3Niv", "ICC" ] <- performance::icc(mn3)[[1]]
summary_stat["MNVar3Niv", "InterClass" ] <- 1 - performance::icc(mn3)[[1]]
summary_stat["MNVar3Niv", "PseudoR2" ] <-
(1 - (summary_stat["MNVar3Niv", "Total_Var" ] / summary_stat["MNNull", "Total_Var" ] ))
summary_stat
ranef <- ranef(mn3)$LIBCOM %>%
mutate(LIBCOM = as.character(rownames(.)),
penteQuality = as.numeric(`Quality`))
mn3_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))
ggplot(data = mn3_ranef) +
geom_sf(aes(fill = penteQuality, colour = Cote)) +
scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
geom_text_repel(data = mn3_ranef[abs(mn3_ranef$penteQuality) > 0.25,],
aes(label = nom_comm, geometry = geometry),
stat = "sf_coordinates",size=3, colour="black")
Résultats différents entre Côte de Beaune (et en particulier Pommard) où les vignobles de grande qualité physique sont globalement associés à des prix plus élevés et Côte de Nuits et Côte d’Or (et en particulier Nuits-St-Georges) où les vignobles de grande qualité physique sont globalement associés à des prix moins élevés…